#include "PhysicPropData.h"

namespace PRS{

	PhysicPropData::PhysicPropData(){
	}

	PhysicPropData::~PhysicPropData(){
	}

	void PhysicPropData::initialize(MeshData *pMData, SimulatorParameters *pSimPar, TMesh* pTMesh_data){
		Swr = pSimPar->Swr();				// Irreducible water saturation
		Sor = pSimPar->Sor();				// Residual oil saturation
		mi_w = pSimPar->waterViscosity();	// water viscosity
		mi_o = pSimPar->oilViscosity();		// oil viscosity
		setInitialSaturation(pTMesh_data,pSimPar);
		initializePressureGradient(pTMesh_data);
		ksModel = pSimPar->ksModel();
	}

	void PhysicPropData::initializePressureGradient( TMesh* pTMesh_data){
		int dim = pTMesh_data->getMeshDim();
		int ndom = pTMesh_data->getNumDomains();
		int k, nrows, ncols;
		for (int dom=0; dom<ndom; dom++){
			pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[dom].getValue(j,0);
				vdata->grad_p = new double[dim];
				vdata->grad_Sw = new double[dim];
				vdata->vel = new double[dim];
				vdata->countedges = 0;
			}
		}
	}

	void PhysicPropData::setInitialSaturation(TMesh* pTMesh_data, SimulatorParameters *simPar){
		bool hasInjectionWell = false;
		int nrows, ncols;
		int ndom = pTMesh_data->getNumDomains();
		// set initial saturation
		for (int i=0; i<ndom; i++){
			pTMesh_data->vertices_list[i].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				VertexData* vdata = pTMesh_data->vertices_list[i].getValue(j,0);
				vdata->Sw = (simPar->isInjectionWell(vdata->flag))?1.0:0.0;
				if ( vdata->Sw > .0 ){
					//printf("Injection well located in node %d Sw = %f flag: %d\n",vdata->ID,vdata->Sw,vdata->flag);
					hasInjectionWell = true;
				}
			}
		}

		// rank with injection well must say to all other ranks that it's OK, a injection well was being informed!
		int well = (hasInjectionWell)?1:0;
		//todo: parallel
		//well = P_getMaxInt(well);
		if (!well){
			throw Exception(__LINE__,__FILE__,"Injection wells are missing!");
		}
	}

	double PhysicPropData::getFractionalFlux(const double &Sw){
		double krw = pow((Sw - Swr)/(1. - Swr - Sor),2);
		double kro = pow((1. - Sw - Swr)/(1. - Swr - Sor),2);
		return krw/( krw + (mi_w/mi_o)*kro );
	}

	double PhysicPropData::getOilFractionalFlux(const double &Sw){
		double krw = pow((Sw - Swr)/(1. - Swr - Sor),2);
		double kro = pow((1. - Sw - Swr)/(1. - Swr - Sor),2);
		return (kro/mi_o)/( krw/mi_w + kro/mi_o );
	}

	double PhysicPropData::get_ksw(const double &Sw){
		switch(ksModel){
		case 1:
			return Sw;
		case 2:
			return pow(Sw,2);
		case 3:
			return pow((Sw - Swr)/(1. - Swr - Sor),2);
		default:
			throw Exception(__LINE__,__FILE__,"Unknown model for relative permeability.\n");
		}
	}

	double PhysicPropData::get_kso(const double &Sw){
		switch(ksModel){
		case 1:
			return 1. - Sw;
		case 2:
			return pow(1. - Sw,2);
		case 3:
			return pow((1. - Sw - Swr)/(1. - Swr - Sor),2);
		default:
			throw Exception(__LINE__,__FILE__,"Unknown model for relative permeability.\n");
		}
	}

	void PhysicPropData::calculateTotalMobility(TMesh* pTMesh_data){
		VertexData *vdata;
		int nrows, ncols;
		int ndom = pTMesh_data->getNumDomains();
		for (int i=0; i<ndom; i++){
			pTMesh_data->vertices_list[i].getSize(nrows,ncols);
			for (int j=0; j<nrows; j++){
				vdata = pTMesh_data->vertices_list[i].getValue(j,0);
				#ifdef ELLIPTIC_PROBLEM
					vdata->lambda_total = 1.0;
				#else
					double krw = get_ksw(vdata->Sw);
					double kro = get_kso(vdata->Sw);
					vdata->lambda_total = krw/mi_w + kro/mi_o;
				#endif
			}
		}
	}
}
